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P5 ■ Abstract 

o , 

We investigate the performance of the GW approximation by comparison to 
exact results for small model systems. The role of the chemical potentials in 
Dyson's equation as well as the consequences of numerical resonance broaden- 
. ing are examined, and we show how a proper treatment can improve computa- 

£NJ | tional implementations of many-body perturbation theory in general. GW and 

exchange-only calculations are performed over a wide range of fractional band 
fillings and correlation strengths. We thus identify the physical situations where 

O- ' these schemes are applicable. 

On 

1 Introduction 

Materials with strong electronic correlation are of considerable interest in solid state 
science, but their computational treatment is notoriously difficult: diagonalizing the 
corresponding Hamiltonians is not feasible for large systems, and the strong correlation 
causes mean-field methods to break down. The collective dynamics of such systems 
can in principle be described exactly by many-body perturbation theory, however. In 
this framework, all exchange and correlation effects are absorbed into the self-energy 
operator S, which may be thought of as a non-local, energy-dependent potential. In 
this paper, we investigate a class of self-energies based on Hedin's GW approximation 
Their diagrammatic representation, which neglects explicit vertex corrections, is 
reminiscent of the Fock exchange potential, but the Coulomb interaction includes dy- 
namic screening. Our first aim is to examine whether numerical improvements can be 
achieved by including further correlation effects in the underlying propagators without 
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changing the diagrammatic form of the self-energy. To assess the performance of these 
schemes, we compare the calculated spectra with exact results for small model systems 
that can still be solved by numerical diagonalization techniques. A similar study was 
recently reported for a two-dimensional Hubbard cluster ||; here we extend that work 
by considering further variants as well as a larger variety of systems, most importantly 
a much wider range of band fillings. Our second aim is to optimize the practical im- 
plementation. To this end, we investigate the treatment of the chemical potentials in 
Dyson's equation and the consequences of resonance broadening in the course of nu- 
merical manipulations. Both points are all too often ignored, but may have significant 
impact on calculated spectra. 

The central quantity of interest is the one-particle Green's function G, whose imagi- 
nary part is directly linked to the spectral function A = ti~ 1 \ Im G\. While many authors 
adopt a momentum and energy representation G(k, u) that follows naturally from the 
band theory of extended systems, we will always consider Grr (u) in real space, which 
is more appropriate for finite clusters. Furthermore, this representation has the advan- 
tage of showing the entire excitation spectrum in each diagonal element of the Green's 
function. While this is not always desirable if one wants to concentrate on the evolution 
of particular quasiparticles, it allows us to judge the performance of any self-energy 
approximation on the basis of a single matrix element. However, we also calculated G 
in reciprocal space for corresponding translationally invariant systems to confirm our 
identification of various spectral features with either quasiparticles or satellites from a 
particular excitation. 

The order of this paper is as follows. Section 2 introduces the model Hamiltonian. 
In section 3 we describe our procedure for obtaining the exact Green's function. In 
section 4 we review the GW approximation and address details of the practical imple- 
mentation. Section 5 lists the approximation variants we consider and gives numerical 
results. Finally, section 6 contains our conclusions. 



The Hubbard model || is the classic example of a Hamiltonian that describes strong, 
short-range electron-electron interaction. It is sufficiently simple to be diagonalized 
exactly for small cluster sizes using standard numerical techniques, yet its physical 
behaviour is non-trivial and reflects many properties of real materials. The model variant 
we employ is a finite chain of M ions with open boundary conditions. Each lattice site 
contains one orbital that can accommodate up to two electrons with opposite spin. 
Doubly occupied orbitals are penalized by a repulsive on-site interaction U, while the 
hopping of transient electrons between neighbouring sites yields an energy gain —t. The 
full Hamiltonian is 



2 Model description 
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Figure 1: Exact Green's function (a) and GW approximation (b) for varying chain length 
M with constant 75% band filling and U = 2. The broad quasiparticle and satellite spectral 
features are insensitive to M, indicating that within the Fermi-liquid limitations of the model 
it is possible to generalize the results reported here. 

where c Rct ,cr ct are the creation and annihilation operators for an electron at site R with 
spin a, nR CT = c Rct cr ct is the particle number operator, and (R, R') indicates a sum over 
nearest neighbours only. We choose the energy norm by setting t = 1. The Hamiltonian 
further contains a local potential Vr that will later serve as a mean-field approximation 
for exchange and correlation. We denote the total electron number by N. 

The properties of the Hubbard model have been thoroughly investigated. In par- 
ticular, the one-dimensional case can be solved analytically using the Bethe ansatz Q 
and, in the limit of infinite chain length, is known to yield a Luttinger-liquid ground 
state. The corresponding Green's function describes a gapless spectrum of bosonic col- 
lective modes involving charge and spin degrees of freedom ||. For finite M, however, 
the renormalized quasiparticle weight factors remain non-zero as long as the Coulomb 
integral U does not exceed a critical value ||, which behaves asymptotically like 1/M. 
In this parameter range the model exhibits Fermi-liquid behaviour. The resolution re- 
quired to differentiate convincingly between Fermi liquids and Luttinger liquids is in 
fact near infinitesimal on the energy scale we consider, and the Lorentzian broadening 
of resonances essentially wipes out features on a genuinely small scale that are of pri- 
mary concern in the distinction between the two. Even so, we have confirmed that all 
systems we study in this paper are comfortably within the Fermi-liquid regime, so that 
the same perturbation methods as for higher dimensions can be applied. 

As we are working with small model Hamiltonians, it is essential to consider the pos- 
sible sensitivity of our results to the parameters in (p. The system size is a particularly 
important aspect. We have calculated the exact Green's function and a correponding 
GW approximation for varying chain length while keeping U = 2 and the fractional 
band filling N/(2M) = 75% constant. In figure [I], as in all later graphs, we show the 
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matrix element ImGi^u) in arbitrary units and align the chemical potentials to facil- 
iate comparison. The number of peaks in the spectral function grows with the chain 
length as expected. However, it is also evident that the qualitative appearance of the 
graphs changes little: the broad quasiparticle and satellite peaks are insensitive to M. 
This becomes even clearer when the integral ImG(u/) do/, which averages over os- 
cillations on a small scale, is considered. We have thus demonstrated that the results 
reported in the following sections are relatively insensitive to the chain length and so 
retain significance beyond the particular model geometry, although of course the for- 
mal extrapolation to M — > oo cannot be made because of the eventual transition to a 
Luttinger liquid. 

3 Exact numerical solution 

The exact one-particle Green's function at zero temperature is defined as 

Grh'(* - = -i(N\T{cnAt)cl a (t')}\N) (2) 

where \N) is the ground state of the interacting iV-electron system, T is Wick's time- 
ordering operator, and CRo-(t) = exp(iHt)cR, CT exp(— iHt) denotes the time-dependent 
wave field operator in the Heisenberg picture. We have suppressed the spin index in G 
because the Green's function is diagonal and degenerate in a. It is convenient to Fourier 
transform (0) to the energy domain and rewrite the Green's function in the form 

GWM = (N\ CRa I c^jN) + (N\cl a 1 c Ra \N). (3) 

oj — + -C/AT oo + rt — rjff 

Here En is the ground-state energy corresponding to \N) and Tl ± denotes the Hamil- 
tonian matrix for N ± 1 electrons. 

The main computational difficulty is that the number of basis vectors of the many- 
body problem grows exponentially with the system size, because it quantifies the (2M!) 
/[(2M — N)\N\] possibilities of distributing iV electrons onto M two-fold spin-degenerate 
orbitals. For a ten-site chain at half filling, the largest model we consider, this implies 
a basis size of 184,756 (although greater chain lengths are feasible if the band filling is 
very small or very large). However, less than 0.01% of the elements of Ti are non-zero, 
so that sparse-matrix techniques may be used to obtain \N) and E^. 

The diagonal elements Grr, which enter the calculation of the electron density and 
other quantities, may be calculated without full matrix inversion by tridiagonalizing 
oj =F i En using the recursion method |7j and starting with the vector c R / CT |iV) 
or c-R a \N). For non-diagonal elements a block recursion must be performed. Neither 
an inversion nor a complete diagonalization would be feasible in terms of computer 
memory: for the example M = N = 10 mentioned above, 210 GB are required to store 
the eigenvectors of 7i ± and 254 GB for the inverse of oo =p 7i ± ± En, although this could 
be somewhat reduced by exploiting symmetry relations. 
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Once the recursion coefficients, i.e., the diagonal elements a n and the off-diagonal el- 
ements b\ of the tridiagonal matrix, are determined iteratively up to a suitable recursion 
depth D, the elements of the Green's function are obtained from 



Even with a basis size of 184,756 about 400 recursions are sufficient, since the number of 
actual spectral features is small compared to the basis size. Ideally a single recursive level 
for each additional peak in the spectrum would be necessary, and indeed we require only 
a few recursions per feature to achieve full convergence in practice. The quasiparticles 
and collective excitations of the iV-electron system are determined by the eigenvalues of 
H.^ and feature as simple poles in the Green's function. For numerical convenience, we 
broaden these sharp resonances into Lorentzians by offsetting the singularities from the 
real energy axis by a distance 5. This procedure does not imply a finite lifetime of the 
excited states. 

4 The GW approximation 

The GW approximation constitutes a diagrammatic expansion of the self-energy that ne- 
glects explicit vertex corrections. However, it includes dynamic screening of the Coulomb 
interaction and is thus capable of describing certain correlation effects. Originally the 
GW self-energy was derived as a first order iterative solution of Hedin's coupled equa- 
tions for the propagators of the interacting many-electron system starting from Hartree 
theory HJ. In this section we briefly review the formalism and subsequently address 
crucial details concerning the computational implementation. 

4.1 The self-energy in the GW approximation 

Starting from a mean-field Hamiltonian that may contain a suitable effective potential 
V, we obtain a zeroth order Green's function 



in terms of the one-particle eigenstates \tp s ) and corresponding energy eigenvalues e s . 
The symbol [Iq denotes the chemical potential. Conventionally, the random phase ap- 
proximation 
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for the irreducible polarization propagator is employed, with t = t — t' and a factor 2 
for spin summation. The dynamically screened interaction is then obtained from 

-l 



W 



RPA, 



00 
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1-p 



RPA, 



u)U 



(7) 



in matrix notation. Finally, the self-energy is given by the expression 



^GW 
J RR' 



T 



to which the GW approximation owes its name, rj denotes a positive infinitesimal. The 
self-energy may be inserted into Dyson's equation to yield the improved Green's function 



G 



GW, 



oo 



1 - G°(u) (V R + £ GW » - V) ' G°(oo 



(9) 



where = U{n^+fijii) indicates the Hartree potential. During the calculation we use 
Fast Fourier Transforms to change between the time and energy domains as appropriate 
in order to avoid costly numerical convolutions. 



4.2 Alignment of the chemical potentials 

Dyson's equation (|9|) combines the equation of motion of the interacting with that of 
the corresponding non-interacting system. The self-energy specifies the deviation of the 
quasiparticle states from the bare electrons and holes upon adiabatic introduction of 
the Coulomb potential. It is thus an equilibrium quantity that should really be calcu- 
lated self-consistently, i.e., the dressed Green's function obtained from Dyson's equation 
is reinserted into the self-energy until convergence has been achieved. This procedure 
is so computationally demanding, however, as to make it unfeasible for large-scale ab 
initio calculations. Furthermore, while the random phase approximation P RPA [G°] by 
construction gives the proper response function of time-dependent Hartree theory, the 
same expression evaluated with a self-consistent Green's function ceases to yield a phys- 
ically meaningful propagator due to the neglect of appropriate vertex corrections. As 
a consequence the spectrum becomes broad and structureless ||. On the other hand, 
the self-energy inherits the chemical potential of the underlying Green's function, and 
a mismatch with that of the dressed propagator may result in wrong time-ordering. 

A possible solution is to use a zeroth order G° to evaluate the self-energy but shift it in 
such a way as to align its chemical potential with that of the dressed Green's function. 
This limited degree of self-consistency, originally suggested by Hedin [l], suffices to 
ensure the correct time-ordering while leaving the response function unchanged and 
thus physically meaningful. From the Fourier transform of we see that shifting G° 
by an amount Co on the energy axis translates into an identical shift of the self-energy 

Eg£(" - a) = f fc RR ,(^ -Co- wOWffiV)^ du/ (10) 

where the contour is closed about the upper half-plane. According to Dyson's equation, 
the chemical potential of the dressed Green's function becomes /i = /io + (V u + J2 GW (fi — 
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Figure 2: The self-consistency shift u, which aligns the chemical potential of the zeroth order 
G° with that of the dressed Green's function derived from it, improves the spectral features 
substantially even for a very weak interaction of U = 1. The exact spectrum is shown for 
comparison. 



Co) — V), where the matrix element is to be taken with the highest occupied quasiparticle 
orbital. In practice we use the corresponding |^ s ) of the non-interacting system. The 
shift is determined so that the chemical potential coincides exactly with that of the 
relocated zeroth order Green's function, hence \i — \i§ +uj. Inserting this relation into 
the previous equation yields the explicit solution 

Q = (y* + _ V y 

Despite its early suggestion, this self-consistency shift is often ignored in ab initio 
band structure calculations, where its impact is normally small. It may substantially 
improve the more sensitive satellite spectrum, however. As an example we consider an 
eight-site chain at half-filling with U = 1. In figure | we display results obtained with 
and without Co. The exact spectrum is shown for comparison. While the shift has little 
effect on the quasiparticle peaks, the improved description of the satellites is evident 
even for such a weak interaction. The results follow our previous demonstration that to 
a large extent Co also restores particle number conservation, which is generally violated 
in non-self-consistent many-body theory ||. All subsequent calculations incorporate Co. 

4.3 Infinitesimal peak broadening 

The GW approximation constitutes a diagrammatic expansion of the true self-energy to 
first order in the screened interaction. The underlying equations <H) to @ may be solved 
in the time as well as in the energy domain. While we employ Fast Fourier Transforms 
to switch between the two as appropriate, most implementations work exclusively in 
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Figure 3: Comparison of spectra calculated with 8 = 0.5 and 8 = 0.02, where the latter was 
broadened afterwards. The smaller Lorentzian width incurs a higher computational cost but 
yields more pronounced spectral features. 

the latter, where all operations are either multiplicative or feasible but computationally 
costly convolutions. It is frequently ignored, however, that the relations are only strictly 
valid in the limit 8^0, because convolutions with a finite displacement of the singular- 
ities from the energy axis will mix the real and imaginary parts of the propagators. As 
a consequence, weak features such as satellites tend to smear out and may even become 
undetectable in the calculated spectrum. As an example, figure ^| contrasts the spectral 
functions from two calculations with 5 = 0.5 and 5 = 0.02, where the resonances in 
the latter were suitably broadened afterwards for the purpose of comparison. We also 
show the exact spectrum. The model specifications are M = 10, N = 14 and medium 
interaction U — 4. The features in the second curve are clearly more pronounced. In 
particular, the reproduction of the satellite at —4 is superior. The downside of using 
small values for 8 is an increased number of sampling points, because the energy resolu- 
tion in a numerical treatment must necessarily exceed the characteristic peak width. In 
the following we always choose 8 as small as computationally possible and only broaden 
the final spectra to achieve a Lorentzian width of 0.5 for visualization. 

5 Comparison of different approximation variants 

In this section we investigate in detail several variants of the GW scheme by direct 
comparison with exact results. The main criterion we apply to judge approximations is 
their ability to reproduce the overall shape of the true spectral function, i.e., the position 
and weight of the quasiparticle excitations as well as their satellites. 
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5.1 Different initial Hamiltonians 

The framework outlined in section 4 allows for considerable freedom in the choice of 
the initial mean-field Hamiltonian. We here consider the following options in order of 
increasing complexity. 

1. Hedin's original iterative derivation suggests that the Hartree approximation V = 
V n with self-consistently determined site occupation numbers should be used, 
although this is rarely done in practice. 

2. The most common choice in ab initio calculations is to start with a self-consistent 
exchange- correlation potential from density-functional (DF) theory, which yields 
the same charge density as the interacting system. We simulate this procedure by 
numerically determining a potential V = V u + V xc that reproduces the occupation 
numbers (nn, CT ) of the exact solution. We still refer to this approach as density- 
functional theory, although it is really a site occupation function(al) theory [ |i~0[ . 



3. We can also evaluate the self-energy (§) with the exact Green's function G, which 
is fully renormalized and contains a background satellite spectrum. While this 
approach itself is of course not directly relevant to practical calculations, it serves 
as an example for implementations that attempt to include a maximum amount 
of many-body features in the initial Green's function. 

It was long regarded as self-evident that including as much information about exchange 
and correlation as possible in the initial zeroth order Hamiltonian will provide an op- 
timal starting point for the iteration. Since the extra computational cost required to 
include a Hartree or local-density mean-field in the Hamiltonian is negligible compared 
to a complete GW calculation, this approach possesses great appeal. As our model is 
numerically solvable, we are in fact able to use the exact Green's function as an extreme 
example of an improved propagator that tries to incorporate as many correlation effects 
as possible up to dynamic renormalization. Besides the conventional random phase ap- 
proximation W rpa [(j°] with density-functional theory as a zeroth order approximation 
we have evaluated the same diagrams using the exact Green's function to obtain the 
more sophisticated H /RPA [G], which contains a rich satellite spectrum. Note, however, 
that this is not the exact screening, since vertex corrections in the polarization are still 
ignored. In figure [| we show the exact spectrum together with the results from the four 
possible combinations of these dielectric functions with G° and G in the self-energy for 
two electrons on a twelve-site chain with U = 2. Evidently the four curves differ very 
little, implying that renormalization does not improve the spectrum, but the computa- 
tional expense is considerably higher if the exact Green's function is used either in the 
screening or the self-energy. 

The rule emerging here is that there is no particular advantage in sophisticated, 
renormalized propagators when one works within the GW scheme. Instead these should 
be consistent with the current level of iteration, because the self-energy (|8]) neglects 
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Figure 4: Using a renormalized Green's function or screening in the GW self-energy without 
vertex corrections fails to improve the spectrum. The numbered curves indicate the combina- 
tions (1) G°W RPA [G°], (2) GW RPA [G°], (3) G°W RPA [G] and (4) GW RPA [G] for two electrons 
on a twelve-site chain with U = 2. 



explicit vertex corrections and so does not become exact when evaluated using the true 
propagators. In general, overrealistic propagators may even cause the approximation to 
deteriorate, because they destroy the balance that exists between the internal diagram- 
matic expansion of the Green's function and screening on the one hand and the vertex 
function on the other [IT]. We thus restrict ourselves to Hartree and density-functional 
theory as zeroth order approximations in the following and present results for these cases 
below. 



5.2 Model dielectric functions 

As the random phase approximation demands inconvenient numerical convolutions in 
the energy domain, practical implementations frequently deviate from the original GW 
scheme by employing alternative model dielectric functions that can be evaluated di- 
rectly ||12|| . In order of increasing complexity, the following options present themselves. 

1. Neglecting screening effects and using the bare Coulomb interaction is particularly 
inexpensive, because no intermediate polarization propagator is required. In this 
case the self-energy 



^RR' 



(12) 



is diagonal and energy-independent. We denote the corresponding Green's func- 
tion, which contains electronic exchange but no dynamic correlation, by G x . If 
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performed self-consistently, this approach is identical to the Hartree-Fock treat- 
ment. 

2. The screening jy RPA [G°] may be calculated in the random phase approximation 
from the zeroth order Green's function in accordance with the original proposition 
of the GW approximation, yielding G GW . 

3. Irrespective of the initial mean- field Hamiltonian, we can employ a realistic screened 
interaction that contains more correlation effects. However, we have already ar- 
gued earlier that a unilateral expansion of the screened interaction will not improve 
the spectrum, and figure ^ gave a numerical example with H^ RPA [G] to this effect. 
Therefore we will not consider this option further. 

In figure ^ we show the exact spectral function for a ten-site chain together with four ap- 
proximations, namely G GW and G x evaluated both with Hartree and density-functional 
theory as zeroth order Hamiltonians. We increase the band filling from 50% in two 
steps to 90% and consider the situation of medium (U = 4) as well as strong {U = 8) 
correlation. Due to the particle-hole symmetry of the Hubbard model, the spectrum 
for band fillings below 50% can be obtained by inflection. For reference, we quote the 
effective potential parameters Vr, for a selected system with 70% band filling and U = 4 
in table [TJ. 

Irrespective of the correlation strength, the occupation numbers are necessarily uni- 
form at half filling due to particle-hole symmetry. Hence in this case the dressed Green's 
functions derived from Hartree and density-functional theory as starting points coincide 
except for a constant shift. The exchange S x is also uniform and causes just another 
constant shift, while the GW approximation improves slightly but still fails to reproduce 
the true spectrum satisfactorily, especially for large U. For instance, the energy gap in 
figure ||(b) is crucially underestimated. On the plus side, the approximations have the 
correct symmetry about the chemical potential. As the band filling increases, the GW 
approximation at first becomes better but again deteriorates for very large band filling. 



The dominant interaction processes in this limit are described by the T-matrix |L3 
which renormalizes the Hartree-Fock potential by including multiple scattering in the 
particle-particle channel to all orders. YP W ignores these diagrams and so does not 
reproduce the spectrum well: it is outperformed even by the unrenormalized exchange 
approximation. The GW approximation works best for intermediate band filling as 
in figures |5](c) and (d), where it yields good results even for strong correlation. In 
this regime there is also little difference between Hartree and density-functional theory 
as starting points, unlike for higher band filling where the GW approximation breaks 
down. The description of quasiparticles is in general superior to that of satellites, which 
are more sensitive to the specific form of many-body interaction processes. The bare 
exchange approximation naturally works best in situations where screening effects are 
negligible, i.e., if the interaction is weak and/or the band filling is very high, such as in 
figure |5|(e). The absence of long-range interaction in the Hubbard model adds to this 
effect. Hartree theory as a zeroth order approximation gives consistently better results 
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Figure 5: The exact Green's function for a ten-site chain compared to the GW approximation 
G GW and the exchange-only G x for different band fillings N/ (2M) and correlation strengths U. 
Where pairs of curves are shown, the solid line refers to a zeroth order density-functional and 
the dotted line to a Hartree Hamiltonian. Some small satellites are marked by arrows. The 
GW scheme performs best for intermediate band filling. The exchange-only scheme yields 
increasingly accurate quasiparticles for high band filling and weak correlation but cannot 
produce satellites. 

Table 1: Effective potential parameters in Hartree and density- functional (DF) theory for 
ten sites with 70% band filling and U = 4. The chain is symmetric about its centre. 



Site index 


1;10 2;9 3;8 4;7 5;6 


Hartree 
DF 


0.26 -0.02 0.06 0.10 0.06 
0.26 -0.33 -0.08 -0.06 -0.15 



Assessment of the GW approximation using Hubbard chains 



13 



than density-functional theory. Being an effective mean field, S x of course cannot pro- 
duce a satellite structure. This is acceptable, however, because collective excitations 
carry little spectral weight in the limits of small U and high fractional band filling. 

6 Conclusions 

We have examined the implementation and performance of the GW approximation by 
comparing exact and approximate spectra for finite Hubbard chains that exhibit Fermi- 
liquid behaviour. The insensitivity with respect to the chain length makes the results 
transferable, although further investigations for higher-dimensional systems will be use- 
ful. Focusing on the computational implementation, we pointed out that the mathemat- 
ical structure of many-body perturbation theory requires a shift aligning the chemical 
potentials in Dyson's equation as well as a minimization of the artificial Lorentzian 
broadening of spectral peaks. In each case model calculations clearly demonstrated 
improvements in the GW approximation, particularly with respect to the satellite spec- 
trum. Next we studied the performance of the GW and the related bare exchange 
approximation for various band fillings and correlation strengths. Our results show that 
the former can yield good results for intermediate band filling even if the correlation is 
strong, while the latter provides a computationally cheaper alternative for weak correla- 
tion and small or high band filling. The description of quasiparticles is generally better 
than that of satellites. We have also shown evidence that renormalized propagators 
will not improve the GW approximation without the inclusion of appropriate vertex 
corrections in the self-energy. 
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